Short-term paleogeographic reorganizations and climate events shaped diversification of North American freshwater gastropods over deep time

What controls species diversity and diversification is one of the major questions in evolutionary biology and paleontology. Previous studies have addressed this issue based on various plant and animal groups, geographic regions, and time intervals. However, as most previous research focused on terrestrial or marine ecosystems, our understanding of the controls on diversification of biota (and particularly invertebrates) in freshwater environments in deep time is still limited. Here, we infer diversification rates of North American freshwater gastropods from the Late Triassic to the Pleistocene and explore potential links between shifts in speciation and extinction and major changes in paleogeography, climate, and biotic interactions. We found that variation in the speciation rate is best explained by changes in continental fragmentation, with rate shifts coinciding with major paleogeographic reorganizations in the Mesozoic, in particular the retreat of the Sundance Sea and subsequent development of the Bighorn wetland and the advance of the Western Interior Seaway. Climatic events in the Cenozoic (Middle Eocene Climate Optimum, Miocene Climate Optimum) variably coincide with shifts in speciation and extinction as well, but no significant long-term association could be detected. Similarly, no influence of diversity dependence was found across the entire time frame of ~ 214 Myr. Our results indicate that short-term climatic events and paleogeographic changes are relevant to the diversification of continental freshwater biota, while long-term trends have limited effect.


Results
The North American fossil freshwater gastropod fauna comprises 606 species in 24 families ranging from the Late Triassic to the Pleistocene (c. 214.0-0.0117 Myr ago). Estimated diversification rates varied little over most of that time interval, especially the extinction rate (Fig. 1). Major peaks in the speciation rate were recorded for the Middle to Late Jurassic (although linked with a wide credible interval) and the Miocene. Elevated speciation rates were also found in the Late Cretaceous and the late Neogene (Figs. 1a, S2a). The extinction rate was elevated in the mid-Cretaceous and reached a minor plateau in the mid-Eocene (Figs. 1b, S2b). The unifactorial models indicated a negative correlation with continental fragmentation, quantified using the Shoreline Development Index (SDI), and a positive correlation with regional temperature (Figs. 2a, 3). There was no significant correlation with diversity. Among all factors, the model with SDI was found to best represent the speciation rate according to Bayes factors (Fig. 2b). Also, the regional temperature model was found to explain variation in the speciation rate slightly better than the constant model (Δ2 log BF = 3.088). For the extinction rate, regional temperature had the lowest median Bayes factors, but the model was not significantly better than a constant model (Δ2 log BF = 0.303; Fig. 2b). Global temperature was positively correlated with trends in the speciation rate but not the extinction rate (Fig. 2a); in both cases, Bayes factors indicated that the models were significantly worse than a constant model (Fig. 2b).

Discussion
Our results indicate that continental fragmentation is a main control for speciation in North American freshwater gastropods over deep time. Regional temperature variation was found to explain speciation rate only marginally better than a constant model, and none of the abiotic or biotic factors showed a significant correlation with extinction rate. However, the lack of any significant correlation is contradicted by a striking coincidence between several shifts and major paleogeographic changes and climatic events (Figs. 1, 3). Below, we discuss the correlation with individual factors, the implications for our understanding of deep-time diversification processes, as well as potential limitations of our interpretations and the use of the selected factors.
Paleogeography: deep-time vs. short-term impact on diversification. Our models suggest that continental fragmentation is a significant control of variation in the speciation rate over long spans of geological time (Fig. 1). This result confirms findings by several other studies suggesting that changes in the paleogeographic setting influence diversification of continental faunas on a large scale 3,21,33,34 . A more fragmented continent leads to poorly connected freshwater ecosystems, which in turn directly affects biogeographic relationships and colonization and speciation processes. Also on a smaller spatiotemporal scale, variation in drainage patterns and hydrological connections plays an important role in the biogeography of freshwater gastropods 35,36 .

Figure 2.
Results of the unifactorial models with speciation and extinction rates. (a) Regional temperature, continental fragmentation (using the Shoreline Development Index) and partly global temperature showed significant correlations with speciation and extinction rates (95% credible interval different from zero). (b) Bayes factors (BF) indicated that the models for SDI and regional temperature are significantly better than a constant model, but only for speciation rate. The plot is truncated at − 100 to focus on the differences in median Bayes factors; see Supplementary Table S1 for the complete data. www.nature.com/scientificreports/ In our case, shifts in the speciation rate coincide strikingly with major changes in paleogeography, more specifically with the opening and closing of intracontinental seaways during the Mesozoic (Fig. 1). The major speciation peak in the Middle to Late Jurassic matches the temporal extent of the Sundance Sea, which covered large parts of today's northwestern United States and southwestern Canada during the Bajocian and Oxfordian stages 31,32 (Fig. 4). The timing of the maximum peak of the speciation rate, however, slightly postdates the end of the Sundance Sea and correlates with the origin of extensive freshwater environments represented in the Morrison Formation (Kimmeridgian-Tithonian) 37 . That formation comprises a diverse set of biofacies, representing different types of freshwater habitats 37,38 , which are summarized under the term "Bighorn wetlands" 37 (here including Lake T' oo' dichi'). A diverse gastropod fauna inhabited this ecosystem, comprising species of Viviparidae, Valvatidae, Lymnaeidae, Planorbidae, Amnicolidae, and Neritidae, amongst others 39 . We hypothesize that the sudden appearance of these new environments provided the ecological opportunities 40 that triggered diversification and resulted in the observed boost in the speciation rate. However, the large credible interval characterizing the peak (Figs. 1a, S2) calls for caution not to overestimate the exact magnitude and duration of the event.
The mid-Cretaceous peak in the extinction rate coincides with the onset of the expansion of the Western Interior Seaway (WIS). This extensive marine incursion intermittently connected the Arctic Ocean with the Gulf of Mexico and the Atlantic Ocean from the late Early Cretaceous (late Albian) to the Paleocene (Danian/Selandian) (Fig. 5) 29,30,[41][42][43] . We hypothesize that the progression of the marine WIS, following a period of extensive freshwater deposition in western North America during the Aptian to earliest Albian 29,30 , caused the extinction of numerous local freshwater faunas. In fact, the loss of this speciose record near the end of the Early Cretaceous is the first in a succession of WIS-influenced faunal turnovers. A second major loss of species occurred at the transition between the Early and Late Cretaceous (late Albian to middle Cenomanian). Most notable for that interval are the Bear River, Cloverly, and Kootenai Formations with faunas dominated by Amnicolidae, Hydrobiidae, Viviparidae, and Pleuroceridae [44][45][46][47] . Following an interval with few faunas (late Cenomanian-Santonian), an abundant continental biota developed in the latest Cretaceous (Campanian and Maastrichtian) with the progradation of freshwater environments 48,49 . The deposits of that time interval (e.g., Foremost, Judith River, St. Mary River, Laramie, Lance, and Hell Creek Formations) are rich in Pleuroceridae, Viviparidae, Physidae, Planorbidae, Hydrobiidae, Ampullariidae, Valvatidae, and other families 47,50,51 . The ongoing faunal turnovers are reflected in the speciation and extinction rates being elevated for millions of years ( Fig. 1).
Despite the coincidence between the peak in extinction and the onset of the WIS, no correlation between extinction rate and continental fragmentation could be detected. This lack of correlation is likely a result of the constant extinction rate over most of the study interval. Generally, we would expect that major paleogeographic changes influence faunal developments due to their impact on the available ecosystems and hydrological connectivity, while smaller or gradual shifts in the continental landscape may be difficult to detect and can cause an overprinting (absent) long-term trend.  53 . Both events triggered shifts in the diversification and/or the migration of mammals, including hominoids 54,55 . Both MCO and MECO were also suggested to be responsible for extinction crises and migrations in European land snails 56 . Conversely, the shifts in speciation and extinction we found in the Mesozoic correlate with general cooling trends (Fig. 3). Despite major climate changes over the study time interval (Fig. 3), prolonged periods are characterized by constant diversification rates (Fig. 1). This is probably the reason why a long-term association explains rate variation only marginally if at all better than a constant model (Fig. 2b), much like for the link between extinction rate and continental fragmentation. Only short-term and rapid climatic events seem to affect diversification of freshwater gastropods, while long-term climate change is apparently of little relevance. This pattern contrasts the recent finding that an interaction of long-term and short-term temperature changes affect speciation rates in marine Phanerozoic biota 57 .
Our result has an immediate implication for the current biodiversity crisis. The current pace of global climate change exceeds by far a natural rate 58 . North American freshwater gastropods are already under severe pressure through habitat destruction or alteration, pollution, and invasive species, amongst other factors 59 . Predicted future shifts in freshwater habitat suitability as a result of climate change across North America will likely increase stress on freshwater biota [60][61][62] . Our results imply that the rapid climate deterioration will come with massive long-term ramifications for diversification in North American freshwater gastropods 63 . www.nature.com/scientificreports/ The second aspect of our results on temperature concerns the finding that global temperature is a poor proxy for deep-time diversification on a continental scale. Although studies have shown a significant relationship between continental diversification trends and global temperature 2 , our analyses indicate that regional estimates should be used instead 52 . Over extended geological time spans, tectonic displacement of continents across different latitudinal/longitudinal ranges involves shifts in regional climatic regimes. This is also the case for the North American continent during the Mesozoic and Cenozoic 27 . Consequently, we advise the use of regional temperature estimates based on raster climate models to assess the impact of temperature variation on diversification in continents or regional settings in general.
Diversity dependence is of little relevance. Diversity dependence seems to have no influence on the diversification of North American freshwater gastropods on an extended temporal scale. There has been a long and intensive debate concerning the influence of diversity dependence and whether there is an upper limit to diversity in general 24,64,65 . Conceptually, diversity dependence refers to the concept that increasing diversity leads to increasing biotic interaction (e.g., competition for limited resources), which in turn causes a decline in the speciation rate and/or a rise in the extinction rate 11,65 . Numerous studies found support for diversity dependence across different species groups, time intervals, and geographic regions 2,4,[15][16][17]65 . Among European fossil freshwater gastropods, diversity dependence was recently reported as the most important control on diversification over short geological timescales 18 . www.nature.com/scientificreports/ The discrepancy among the results found by previous studies, also concerning the question if diversity dependence or upper limits to diversity exist at all, is likely owed to the different temporal, spatial, and taxonomic scales applied 65 . Our results indicate that diversity dependence is only of limited relevance to freshwater biota on a continental scale over extended geological time spans. In contrast, biotic factors are presumably more important at local and regional scales (e.g., in single drainage basins or biogeographic regions) and over shorter time spans 18,23 . The absence of diversity dependence of the extinction rate suggests that North American freshwater gastropod faunas were far from approaching equilibrium diversity.
Limitations. The uneven sampling (Figs. S1, S3) and variable preservation over geological time and across the North American continent complicate reconstructions. Fossil data are derived primarily from the western plains and intramontane basins of the United States and Canada, while other parts of North America are rather poor in freshwater deposits and freshwater mollusks 66,67 (Fig. S1). However, since the links to environmental changes concern globally relevant climatic events (MECO, MCO) and paleogeographic reorganizations concentrated in regions with rich faunas, we do not expect a strong bias in our reconstructions. Nonetheless, potential environmental events that may influence diversification in regions where no or few faunas are preserved are not covered by this study. This concerns, for example, large parts of eastern North America or the regions where marine deposition prevailed over prolonged periods (e.g., those covered by the Western Interior Seaway), which left limited space for freshwater ecosystems and their faunas to unfold. Similarly, certain time intervals lack data, such as the period covered by the Sundance Sea or parts of the Early Cretaceous (Fig. S3a). This results in wider credible intervals in speciation and extinction as well as preservation rates (Figs. 1, S2, S3b), which is why associated peaks or apparently constant rates should be interpreted with caution. Future investigations should attempt to incorporate spatial biases in the analyses [68][69][70] .
Finally, it is important to note that, in contrast to Europe, North America harbored only few long-lived lakes, in which species diversity and endemism tends to be high. Long-lived lakes typically produce extensive sedimentary deposits that have a higher chance of preserving mollusk faunas than their short-lived equivalents or rivers 71 . A prominent example is Late Jurassic Lake T' oo' dichi' , an alkaline-saline wetland-lake system that existed for approximately 2 Myr in the southern part of the Bighorn wetlands and deposited parts of the Morrison Formation 37,72 . A Cenozoic representative is the lake system in the Uinta-Piceance-Greater Green River Basin system, which contains lacustrine deposits that span several million years in the early Eocene 73 . Instead, North America's fossil freshwater faunas are more commonly found in near-shore and riverine settings as well as short-lived lacustrine systems. Also, the poor preservation of shells in the many North American freshwater systems limits the actual known diversity 74,75 . The difference between North American and European freshwater faunas raises the interesting question whether the absence of long-lived lakes in North America could have depressed both speciation and extinction rates relative to the pattern in Europe, where the origin and demise of large ancient lakes greatly increased speciation and extinction 33,76,77 .

Conclusion
Our results highlight the relevance of abiotic factors for diversification of continental freshwater biota and the timescales at which they act. Diversification rates of North American freshwater gastropod faunas are constant over long periods of time. Continental fragmentation-especially the changes in hydrological connectivity of freshwater ecosystems it entails-and regional temperature showed a significant long-term effect on the speciation rate, while variation in the extinction rate could not be explained by any factor. Diversity dependence showed no correlation at all. The apparent long-term influence of the abiotic factors, however, seems only a combined result of the impact of abiotic factors during selected periods or even short-term events. Notable peaks in the speciation and extinction rates coincide with major paleogeographic reorganizations in the Mesozoic, i.e., the formation of the Bighorn wetlands following the Sundance Sea (Late Jurassic) and the Western Interior Seaway (late Early Cretaceous to Paleocene), as well as with climatic events in the Cenozoic, i.e., the Middle Eocene Climate Optimum and the Miocene Climate Optimum. Especially the formation and retreat of the Mesozoic seaways coincided with extensive adjacent freshwater deposition, which likely provided impetus for the observed shifts in diversification. Considering the episodic nature of the climatic and paleogeographic changes coinciding with rate shifts, compared to long intervals without detectable shifts in diversification, statistical analyses based on long-term trends may fail to detect short-term influence. As such, our working hypothesis is only partly accepted-abiotic factors, and especially continental/hydrological fragmentation, are more important for diversification processes than biotic factors, but a significant long-term effect could not be supported. Future research may lead to the identification of specific time intervals or geographic regions to disentangle the diversification processes at different spatiotemporal scales.
We compiled the fossil species occurrences of Late Triassic to Pleistocene North American freshwater gastropods from the literature 18 (Fig. S1). For systematic classification (genus and family), rank (species or subspecies) www.nature.com/scientificreports/ and taxonomic status (valid or synonym), we followed the latest opinion as published in the literature. Likewise, the stratigraphic attribution of the fossil-bearing strata was taken from the most recent geological assessments to avoid a bias of the results from outdated chronostratigraphy. Species with uncertain identification, nomina dubia, nomina nuda, and taxa inquirenda were excluded from the analyses, while "good" taxa with invalid names (e.g., being junior homonyms pending revision) were included. Contrasting other studies on long-term diversity or diversification trends that focus on the genus level 99-101 , we use the species as the basic unit for the analyses for two reasons. First, using genera in evolutionary studies has been shown to severely bias the results in some cases 102 . Second, the genus classification is difficult for many fossil freshwater gastropod species given the sometimes poor preservation (especially in Mesozoic strata) and the paucity of diagnostic characters. This typically leads to classification in overly broad (and often extant) genera that results in unrealistically long stratigraphic ranges.
The geographic delimitation of North America through geological time follows the present-day borders of countries and ranges from Alaska, Canada, and Greenland in the North to the Isthmus of Panama in the South. Excluded are Iceland and Svalbard, which are presently associated with Europe, but were connected to the North American continent at times 27 . Despite major paleogeographic changes affecting the connectivity of North America and adjacent continents since the Triassic 27 , assemblages known from coeval deposits of other continents bear little resemblance to the studied faunas 103 . Except for very few widespread species in the Pleistocene, which have a (nearly) Holarctic distribution [e.g., Aplexa hypnorum (Linnaeus, 1758) 104,105 ], the North American fossil freshwater gastropod fauna is distinct. Alleged records of species originally described from North America in Asian fossil deposits 106  Diversification rate analyses. Diversification rates were inferred with the free software PyRate v. 3 108 (https:// github. com/ dsilv estro/ PyRate) following the protocol of Ref. 109 . To take age uncertainties into account, we ran the analyses for 100 randomized data sets, resampling the age of each species occurrence evenly within its fossil range. The algorithm jointly models fossil species occurrences as a result of speciation, extinction, and preservation, while accounting for sampling heterogeneity through geological time. Sampling is modeled as a time-variable Poisson process according to stratigraphic epochs. The Reversible Jump Markov Chain Monte Carlo (RJ-MCMC) algorithm was used to infer time-varying speciation and extinction rates, including the estimation of the number and temporal placement of rate shifts over time.
The basic analysis setting involved 50 million MCMC generations with a sampling frequency of 5000. Output files yielding a low effective sampling size (ESS < 100) were re-analyzed with 100 million MCMC generations but at the same sampling frequency to achieve convergence. ESSs were calculated with the package 'coda' v. 0.19-3 110 for the statistical programming environment R v. 4.1.2 111 . In all cases, the initial 20% of the samples were discarded as burn-in. Rate plots were created at a temporal resolution of 0.1 Myr and are based on 100 randomly sampled generations per replicate.

Potential controls on diversification.
To test for a link between diversification rates and paleogeography, climate, and species diversity, we employed the univariate birth-death model implemented in the PyRate package 112 . As a measure for continental fragmentation we used the Shoreline Development Index (SDI). This index was originally introduced to quantify the outline complexity of lakes and describes the perimeter (P) deviation from a perfect circle with the same area (A) as the lake 113 : We reconstructed the area and perimeter of the North American continent through time based on latest paleogeographic reconstructions 27 . The study time interval (Late Triassic-Pleistocene) includes the opening of the northern Atlantic Ocean, when North America adjoined to large parts of South America, Africa, and Europe during the Late Triassic and Jurassic. Where seaways were absent, the border of North America was defined based on boundaries of tectonic plates and reconstructed country borders, as employed by the PALEOMAP project 114 . Boundary polygons were created in GPlates v. 2.2 115 . Using the R packages 'raster' v. 3.4-5 116 , 'rgeos' v. 0.5-5 117 , and 'sf ' v. 0.9-7 118 , paleogeographic maps were cropped based on the boundary polygons and converted themselves to polygons. These were subsequently simplified using the function 'gSimplify' of the R package 'rgeos' with a numerical tolerance of 10 and afterwards smoothed using the package 'smoothr' v. 0.2.2 119 and a Gaussian kernel regression with a smoothness factor of 3. Islands and holes in polygons greater than 5000 km 2 were deleted subsequently with 'smoothr' . The final polygons defining the North American continent at a given time slice were unified, and areas and perimeters were calculated with the package 'geosphere' v. 1.5-14 120 based on a WGS84 ellipsoid (Fig. 3a,b).
The simplification and smoothing procedure hardly affects the area approximation but is particularly relevant for the reconstruction of comparable perimeters through time as past outlines are less certain than present-day ones. However, despite these efforts, the polygon for today's North America is still much more accurate and thus has a higher outline complexity. To avoid any bias of the analyses from these differences in resolution, we omitted the recent time bin and ran the analyses (for all abiotic and biotic factors) for the time interval 5-215 Myr ago. In doing so, we also limited a potential pull-of-the-recent effect.
For temperature, we used both global and regional estimates to test for a variable influence of the spatial scale of observation. We focus here on regional temperature, considering that climate variation across the North American continent is more likely to directly influence diversification than global temperature 52  www.nature.com/scientificreports/ also include global temperature as it is a commonly used measure in diversification studies. Using both factors allows testing for differences in the impact of regional and global climate. Global average air temperature derived from the latest model of Ref. 28 on a 1-Myr resolution (Fig. 3d). Regional temperature was reconstructed based on global temperature raster data on a 1° × 1° spatial resolution. These rasters derive from the PALEOMAP Project, which has produced paleogeographic maps at 5-Myr intervals 121 and has assembled related temperature data based on the HadleyCM3 paleoclimate simulations 122 , which use estimates of the past concentration of atmospheric CO 2 123 . It is a well-known observation that global climate models tend to produce cooler temperatures at higher latitudes, especially for hot house time intervals. In order to correct this bias, we adjusted the latitudinal temperature gradient using the Phanerozoic climate model 28 . This resulted in an average increase in temperature of c. + 5 °C at 60° latitude. North American average ("regional") temperatures were calculated based on the intersection of the global temperature raster data and paleogeographic reconstructions 27 with the R package 'raster' (Fig. 3d). The data for SDI, global, and regional temperature were centered and scaled by their standard deviations to facilitate comparison and interpolated to 0.1 Myr increments.
To account for age uncertainties of the fossil record, estimated times of origination and extinction for each species were extracted once for each of the 100 replicates (see previous chapter). The exponential correlation model was used for all three predictor variables and compared to a constant model, where rates are treated as invariable through time 112 . The unifactorial diversification models were ranked by their fit as quantified through marginal likelihoods K, which were inferred through thermodynamic integration in PyRate 108 . This method approximates the marginal likelihood of the data via ten MCMCs with increasingly altered acceptance ratio 124 , each running for 1,250,000 generations with a sampling frequency of 5000 and a burn-in of 10%. Output files yielding an ESS < 100 were re-analyzed with 2,500,000 or 3,000,000 MCMC generations but at the same sampling frequency to achieve convergence. Median Bayes factors (BF) comparing marginal likelihoods of alternative models across the 100 replicates were evaluated, whereas Δ2 log BF > 2 was considered as evidence for significant differences 125 . Correlations with biotic and abiotic factors were found significant if the 95% credible interval did not include zero 108 .

Data availability
The data that support the findings of this study are available in the supplementary material of this article. The dataset is available at https:// doi. org/ 10. 22029/ jlupub-466. www.nature.com/scientificreports/